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Abstract. We explore a new type of discretizations of lattice dynamical mod¬ 
els of the Klein-Gordon type relevant to the existence and long-term mobility 
of nonlinear waves. The discretization is based on non-holonomic constraints 
and is shown to retrieve the “proper” continuum limit of the model. Such dis¬ 
cretizations are useful in exactly preserving a discrete analogue of the momen¬ 
tum. It is also shown that for generic initial data, the momentum and energy 
conservation laws cannot be achieved concurrently. Finally, direct numerical 
simulations illustrate that our models yield considerably higher mobility of 
strongly nonlinear solutions than the well-known “standard” discretizations, 
even in the case of highly discrete systems when the coupling between the ad¬ 
jacent nodes is weak. Thus, our approach is better suited for cases where an 
accurate description of mobility for nonlinear traveling waves is important. 


1. Introduction. Discrete solitons such as kinks and breathers exist in many phys¬ 
ical systems, including Josephson juctions, electrical circuits, and granular crystals 
[1, 2, 3]. In the models that describe these systems, the mobility of soliton solu¬ 
tions, and specifically kinks - which will be the focus of the present work-, is an 
important and well-studied property [4, 5]. For instance, kinks in extemely dis¬ 
crete systems need to overcome the celebrated Peierls-Nabarro barrier, while they 
undergo energy radiation which leads to their trapping in the lattice [9]. More re¬ 
cent work revealed the existence of kinks that do not experience the action of the 
Peierls-Nabarro potential [10]. 

For a class of Klein-Gordon discretizations [5] it has been shown that energy and 
linear momentum cannot be both conserved. This is problematic in cases when the 
objective is the discretization of the continuum version of equation, which conserves 
both momentum and energy. Going beyond the properties of isolated kinks, the 

2010 Mathematics Subject Classification. Primary: 37K60; Secondary: 81Q05. 

Key words and phrases. Klein-Gordon models, discrete solitons, kinks, Peierls-Nabarro barrier, 
conservation laws. 


1 



2 


P.G. KEVREKIDIS AND V. PUTKARADZE AND Z. RAPTI 


above result has consequences for the studies of soliton collisions in Klein-Gordon 
lattices. In general, the mobility of kinks in discrete manifestations of the Klein- 
Gordon equation may vary from that in the continuous version as well as among 
different discretizations. 

In this work, we present a class of discretizations based on a non-holonomic con¬ 
straint that preserve the momentum and alleviate the problem of decreased mobility 
in kink solutions. Specifically, as we enforce the conservation of momentum it in 
turn implies an affine constraint in kink velocities, which is of the non-holonomic 
type. The extra term added in the discretization due to the non-holonomic con¬ 
straint, vanishes in the continuum limit, hence demonstrating that the discretized 
version of the equation converges to its continuous counterpart. Our numerical re¬ 
sults validate our theoretical considerations and we present multiple instances where 
the additional constraint that we impose increases the mobility of the kinks. 

We shall note that for generalized non-holonomic constraints that are not linear 
in velocities, one has to be careful in applying the variational principle. This was 
demonstrated in, for example, [6, 7] for the case of stabilizing servomechanisms, 
where an alternative formulation through Poisson structures has been derived. The 
discrete equivalent of a conservation law for continuous system (1) below, enforced 
by the constraint, cannot be readily interpreted as a consequence of a continuous 
symmetry through a Noether theorem. Thus, applications of alternative methods 
of non-holonomic mechanics [8], in addition to the Lagrange-d’Alembert’s principle 
developed here, are of particular interest and will be considered in future work. 

This article is structured as follows. In section 2 we describe the Klein-Gordon 
model that will be studied and present the theoretical framework of our method. In 
section 3 we present numerical results that corroborate our theoretical predictions. 
Finally, in section 4 we discuss our findings and offer insight for future directions. 

2. Model Theoretical Setup and Proposed Discretizations. The models 
that we will consider herein will involve discretizations of continuum partial dif¬ 
ferential equations of the nonlinear wave type: 

~ ’^XX ^ (^)' ( 1 ) 

There are numerous key model examples of this class of equations, such as the 
sine-Gordon [1] or the (j)^ model arising from classical field theory and high-energy 
physics [2], as well as variants thereof. 

The corresponding discrete models are of considerable interest in their own right 
as models of coupled torsion pendula, dislocation dynamics, Josephson junctions, 
electrical circuits, among numerous other settings [1, 2, 11, 12]. The generic discrete 
analogue of these models is derived by using a standard lattice discetization of the 
second space derivative as follows: 

iln — -\- Uji — l V (Un) (2) 

The field is defined at the integer nodes n and e denotes the coupling constant 
between adjacent sites, which could be interpreted as e = Ijdx^, with dx being the 
spacing between the adjacent sites (see below). 

Note that lattice model (2) is Hamiltonian with 

H = ^ ^ 2 ^””^ 2 ~ ^{Un) = 'y ( T~Lm 

n n 

where T-Ln is the energy density. 


( 3 ) 
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In that context, an alternative limit, namely e = 0, the so-called anti-continuum 
limit becomes particularly important and provides a valuable tool for assessing 
the potential stationary or time-periodic states of the model and their respective 
spectral stability properties [13, 14, 15]. 

One of the “pathologies” of such discrete models is that as the coupling param¬ 
eter e decreases, the mobility of the nonlinear waves accordingly weakens [4]. The 
decreased mobility stems, to a considerable degree, from the fact that the momen¬ 
tum P = J utUxdx is no longer a conservation law in the discrete case of Eq. (2), 
even though it is in the corresponding continuum limit of e —>■ oo. This curious fact 
has spurred a considerable volume of recent work (for a partial summary, see [16]) 
focused on providing alternative discretizations that would enable higher mobility of 
the waves, even at small e. Our aim in the present work is to propose a novel class 
of such discretizations, based upon the imposition of a non-holonomic constraint 
providing an exact conservation of the discrete momentum. 

Indeed, the continuum equation (1) provides a conservation of the momentum 
quantity Pc = f utUxdx. On the other hand, after defining the analogue of discrete 
momentum for (2) as 

P = '^U„(U„+1 - Un-l) (4) 

a quick check shows that such discrete momentum is not conserved [17]. The goal of 
this paper is to enforce the conservation of (4) using the methods of non-holonomic 
mechanics. The equation (2) can be obtained as the Euler-Lagrange equation with 
the Lagrangian 

through the critical action principle S f Ldt = 0 on arbitrary variations 6un (t) with 
appropriately vanishing boundary conditions at the end time points. Next, enforcing 
the conservation law (4) leads to an affine constraint in velocities of the type 

an{u)un = b{u) with a„(u) = Un+i - Un-i , b{u) = P (6) 

n 

using the notation u = (ui, M 2 , ... , Un)- Since this constraint involves velocities in 
a non-trivial way, i.e., that constraint cannot be integrated to contain coordinates 
only, the constraint is non-holonomic, and linear in velocities. The usual way of 
dealing with this type of constraints is through the Lagrange-d’Alembert’s principle 
[18]. For the constraints of the type (6), the variations Sun are no longer arbitrary, 
but must satisfy the linear relations 

'^an{u)SUn='^iUn+l-Un-l)5Un=0. (7) 


Notice that we do not enforce (6) directly, as this would correspond to the so- 
called vakonomic approach, relevant for e.g. the control theory but normally not 
leading to correct expressions for mechanical systems, as described in [19]. On 
the contrary, Lagrange-d’Alembert’s approach (7) was consistently shown to be 
physically appropriate. For more in-depth discussion of various approaches to non- 
holonomic systems, we refer the reader to [18, 19]. 

If r is the Largange multiplier enforcing the constraint (7), the critical action 
principle coupled with the constraint gives 



d dL 
dt dUn 


e {Un+l + Un-1 - 2Un) + V {Un) - T (u„+i - Un-l) 


SUn = 0 , 
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provided Sun vanish at t = to,i so no boundary integration terms appear in the 
equations. Since Sun is arbitrary, the equations of motion conserving momentum 
are 


Un= e {Un+l + Un-l - 2Un) - V\Un) + r{Un+l - U„_i) . 


( 8 ) 


Note the new term proportional to r and enforcing the constraint on the right hand 
side of (8), as compared to (2). The physical meaning of r is the (generalized) 
constraint force. In order to find r, we compute the derivative dP/dt = 0 and set 
it equal to zero. We find that while the terms associated with e vanish telescopi¬ 
cally [17], the nonlinear term does not. This leads to the explicit condition for r 
as 

^ ^ Y.nVM{Un+l-Un-l) 

X/n ~ Un — l) 

Making this selection for r, multiplying both sides of Eq. (8) by Un, and summing, 
we observe that the derivative of the (former) Hamiltonian of Eq. (3) is not zero 
anymore, but is given by the following condition: 


dH 

dt 


= rP. 


( 10 ) 


Importantly, this identity suggests that for generic initial data, such that P ^ 0, 
the Hamiltonian will not be conserved. This is due to the constraint (6) being 
non-homogeneous in velocities for P ^ 0. It is relevant to point out that such an 
“incompatibility” of the momentum and energy conservation laws for other broad 
classes of discretizations of the Klein-Gordon dynamical lattices have also been 
reported previously [5]. 

An additional comment about this non-holonomic discretization is that the extra 
term introduced by the non-holonomic constraint vanishes at the continuum limit 
e —>■ oo, ensuring that Eq. (8) converges to the proper continuum limit of Eq. (1). To 
see this, we consider the continuum analogue of the new r{un+i —Un-i) term, which 
is Ux(f V'{u)uxdx/ f u^dx). The numerator of the term in the parenthesis can be 
integrated to yield V{u 2 ) — V{ui) where U 2 and ui are, respectively the asymptotic 
values of the field at ±oo. Assuming that we examine the dynamics of homoclinic 
orbit or of a heteroclinic orbit in a symmetric double well (as is typically the case 
in the examples of interest mentioned above), this quantity vanishes identically and 
the discretization reverts to its proper continuum limit. We will return to this point 
and corroborate it through our numerical computations below. A final remark about 
the nature of our discretization procedure, different from all other schemes that we 
are aware of, is that our method is inherently nonlocal since the computation of r 
in (9) entails a summation over the entire lattice. 


3. Numerical Computations. To provide a comparison between the properties 
of the different discretizations, we examine the dynamics of Eq. (2) vs. that of 
Eq. (8). We use an initial condition in the form of a prototypical nonlinear excitation 
such as the kink. In the continuum limit, the latter stationary solution (which 
can be boosted to any sub-luminal speed via Lorentz transformation) is of the 
form u{x) = tanh(x). This is the case for the (j)'^ model with a potential V{u) = 
{l/2){v? — 1)^. For the discrete system, we have used the initial condition it„(0) = 
tanh(a:) = tanh(n x dx), where dx is used to signal the spatial discretization step 
(see also below), approximating the continuous solution on a discrete set of points. 
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Figure 1. Evolutionary dynamics of the “standard” discrete 
model of Eq. (2) with dx = 1/4, hence e = l/dx'^ = 16 for 
V{u) = — 1)^ i.e., the model. The top left panel repre¬ 

sents the space-time (n — t) contour plot of the energy density (see 
Eq. (3)), while the top right represents the corresponding contour of 
the field u„(t). The bottom left illustrates the “near-conservation” 
for this (approaching the continuum) case for the linear momentum 
of Eq. (4), while the bottom right illustrates that the model is en¬ 
ergy preserving up to fluctuations of size of 10“^"^ due to numerical 
roundoff errors. 


Eurthermore, Eqs . (2) and (8) are integrated using a fourth-order Runge-Kutta 
scheme with a temporal discretization step of the typical size dt = 0.00025. 

In addition, we enforce an initial momentum of Pfix = 0.5 via the definition 
of Eq. (4). The simplest way that we envisioned for achieving this was to assume 
all sites as having vanishing velocity except for the central one of the chain, say 
no, and then using = Pfix/(uno+i — Uno-i)- This prescription provides the 
given amount of momentum to the central portion of the kink and we expect that 
it should set the kink in motion. 

This is exactly what we observe in Eig. 1, which is close to the continuum limit 
of e = 1/dx^ —>■ oo. Here, as usual dx is effectively the corresponding spacing 
between the lattice sites, leading as a finite-difference three-point approximation to 
the second derivative. In particular, in this figure dx = 1/4 or e = 16 is used. The 
top left panel illustrates the traveling kink of the standard discretization of Eq. (2). 
The proximity of the model to its continuum analogue coupled with the imparting 
of momentum via the initial condition leads the kink to travel through the chain 
with an apparently constant speed (evident in the constant slope of the two contour 
plots of the top panel, the left one illustrating the energy density of Eq. (3) and the 
right one the actual field Un{t)). Given that the model of Eq. (2) is Hamiltonian, our 
explicit fourth-order time integration (Runge-Kutta) scheme conserves the initial 
energy (of 0(1)) up to a few parts in 10“^^ verifying its conservation for all practical 
purposes. On the other hand, the momentum fluctuates in the third decimal digit 
indicating its non-preservation by the scheme; instead, it is only approximately 
conserved because of the proximity to the continuum limit where the continuum 
analogue of Eq. (4) can be shown to be preserved by Eq. (1). 

On the other hand, we can see that the fluctuations in this quantity {P) are 
far more dramatic in the considerably more discrete case e = 1, corresponding to 
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Figure 2. Same as the previous plot, only now we consider a 
highly discrete case dx = e = 1. Here, we observe the detrimen¬ 
tal effect of the energy-conserving discretization with regards to the 
mobility of the kink initial condition. The strong discreteness/weak 
coupling precludes the kink from moving through the lattice and 
the initial momentum (which as can be seen is far from being con¬ 
served now) is expended in a jerky motion with a range of a few 
discrete sites (trapping of the kink). 


dx = 1, shown on Fig. 2. There, the fluctuations in P are of 0(1) while the energy 
is still conserved with the same precision as before. In addition, the motion of the 
kink is fundamentally different, featuring a jerky back and forth trajectory between 
a few sites around the origin. Thus, in this case, the kink appears trapped due to 
the weak inter-site coupling which does not allow the momentum to translate itself 
into mobility as it would in the corresponding continuum limit. 

We now turn to the corresponding dynamics of the newly proposed model of the 
form of Eq. (8). The case of dx = 1/4 is again shown in Fig. 3. The fundamental 
difference in this case is that the momentum P of Eq. (4) in the bottom left panel 
is conserved to a similar accuracy as the energy in the previous model, as expected 
by our construction. On the contrary, the energy of the bottom right panel varies 
at the third significant digit remaining roughly constant solely due to the proximity 
to the continuum limit, where the analogue of Eq. (3) is conserved for Eq. (1). 

The main advantage of our discretization method is evident in Eig. 4. More 
specifically, as the example oi dx = e = 1 shown in Fig. 3 illustrates, the discretiza¬ 
tion of Eq. (8) presented here can accurately describe traveling waves in the case 
of strong coupling. Fig. 3 shows that the momentum preserved up to numerical 
accuracy, as expected, while the energy fluctuations become more significant, re¬ 
vealing the absence of energy conservation already in the second significant digit. It 
should also be mentioned that we have also separately verified Eq. (10) numerically, 
although we do not present it here. Hence, our principal numerical conclusion is 
that the models of Eq. (8) ensure a higher mobility of their corresponding coherent 
structures compared with the standard discretization of Eq. (2). This comes at the 
price of the incorporation of the nonlocal effect through the summation that ensures 
the conservation of the linear momentum. 

Lastly, we have confirmed that the relevant discretizations are indeed properly 
converging to the continuum limit as da; 0. A typical example of this approach 
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Figure 3. Same as Fig. 1 but now for the momentum conserv¬ 
ing discretization of Eq. (8). This is the near continuum case of 
dx = 0.25 which sustains motion of the kink with constant speed. 
The momentum of the bottom left is conserved up to numerical 
accuracy, while the energy of Eq. (3) in the bottom right is only 
approximately conserved, due to the proximity to the energy con¬ 
serving continuum limit. 


is shown in Fig. 5. In particular, the figure showcases the evolution of the quantity 
r{t), which constitutes the prefactor of the newly proposed nonlocal term arising 
through the non-holonomic constraint. It can be seen that as we go from the weakly- 
coupled strongly discrete case of the left panel (for dx = 1) to the right panel, closer 
to continuum case of dx = 0.25, the value of r(t) decreases by more than an order 
of magnitude 


4. Conclusions and Future Challenges. In the present work, we have revisited 
the theme of discretizations of continuum nonlinear partial differential equations 
of the wave type. We have focused on the example of Klein-Gordon equations and 
offered a novel class of schemes which can be summarized as follows: (a) they can be 
straightforwardly designed to conserve linear momentum; (b) when they do so, they 
do not in general conserve the energy, (c) they revert to the continuum PDE in the 
appropriate limit of the lattice spacing tending to zero; and (d) they are based on 
the imposition of a non-holonomic constraint and thus are fundamentally different 
than previously proposed schemes due to their nonlocal nature. We have argued 
that these models have intriguing characteristics, such as the ability to impart 
high mobility to localized coherent states, a feature corroborated via numerical 
computations in the case of kinks of the (f>‘^ model. This feature is preserved even 
in the regime of high discreteness/weak coupling where the standard discretizations 
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Figure 4. Same as Fig. 4 but now for e = 1 (corresponding to 
dx = 1). Notice that despite the high discreteness of the model 
and the weakness of the coupling, the mobility of the state is well 
preserved. 



( 




Figure 5. Corroboration of the convergence of the discrete model 
of Eq. (8) to the continuum limit of Eq. (1). The time evolution 
of the quantity r(t) is shown, in the left panel for dx = 1 and in 
the right for dx = 0.25. It can be seen that the quantity assumes 
significantly smaller values in the latter case in comparison to the 
former, and indeed this trend continues for decreasing values of dx, 
with r —0, as dx —>■ 0. 


under the same initial data produce kinks that are trapped within the Peierls- 
Nabarro energy barrier and cannot move through the chain. 

This study opens an interesting new direction for the investigation of the rele¬ 
vance of such schemes from a numerical perspective and also for a more detailed 
understanding of their properties. While in the present initial communication we 
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have focused on the dynamics of kinks, it would be relevant to explore correspond¬ 
ing structural features for discrete breather type solutions [20] . Furthermore, while 
in the present work we have focused on the one-dimensional domains only, it would 
be particularly interesting to consider generalizing such models in higher dimen¬ 
sions. The latter have progressively become a theme of considerable attention over 
the past few years [21, 22] and hence it would again be particularly relevant to 
compare the standard higher dimensional discretizations, with the ones produced 
by the multi-dimensional generalization of the non-holonomic constraint concepts 
presented herein. Such studies are presently in progress and will be reported in 
future publications. 
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